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Abstract 

The Antarctic Muon And Neutrino Detector Array (AMANDA) is a high-energy 
neutrino telescope operating at the geographic South Pole. It is a lattice of photo- 
multiplier tubes buried deep in the polar ice between 1500 m and 2000 m. The 
primary goal of this detector is to discover astrophysical sources of high energy neu- 
trinos. A high-energy muon neutrino coming through the earth from the Northern 
Hemisphere can be identified by the secondary muon moving upward through the 
detector. 

The muon tracks are reconstructed with a maximum likelihood method. It models 
the arrival times and amplitudes of Cherenkov photons registered by the photo- 
multipliers. This paper describes the different methods of reconstruction, which 
have been successfully implemented within AMANDA. Strategies for optimizing the 
reconstruction performance and rejecting background are presented. For a typical 
analysis procedure the direction of tracks are reconstructed with about 2° accuracy. 

Key words: AMANDA, track reconstruction, neutrino telescope, neutrino 
astrophysics 

PACS: 95.55.Vj, 95.75.Pq, 29.40.Ka, 29.85,+c 



1 Introduction 

The Antarctic Muon And Neutrino Detector Array [1], AMANDA, is a large volume 
neutrino detector at the geographic South Pole. It is a lattice of photo-multiplier tubes 
(PMTs) buried deep in the optically transparent polar ice. The primary goal of this 
detector is to detect high-energy neutrinos from astrophysical sources, and determine 
their arrival time, direction and energy. When a high-energy neutrino interacts in the 
polar ice via a charged current reaction with a nucleon N: 

v t + N -> £ + X , (1) 

it creates a hadronic cascade, X, and a lepton, £ = e,/i,T. These particles generate Che- 
renkov photons, which are detected by the PMTs. Each lepton flavor generates a different 



2 



signal in the detector. The two basic detection modes are sketched in figure 1. 




Figure 1. Detection modes of the AMANDA detector: Left: muon tracks induced by 
muon-neutrinos; Right: Cascades from electron- or tau- neutrinos. 

A high-energy charged current interaction creates a muon, which is nearly collinear 
with the neutrino direction; having a mean deviation angle of ip = 0.7° x (E u /TeV)~ ' 7 [2], 
which implies an accuracy requirement of < 1° for reconstructing the muon direction. 

The high-energy muon emits a cone of Cherenkov light at a fixed angle 9 C . It is determined 
by cos# c = (n ■ f3)~ 1 , where n ~ 1.32 is the index of refraction in the ice. For relativistic 
particles, (3 ~ 1, and 9 C ~ 41°. The direction of the muon is reconstructed from the time 
and amplitude information of the PMTs illuminated by the Cherenkov cone. 

Radiative energy loss processes generate secondary charged particles along the muon 
trajectory, which also produce Cherenkov radiation. These additional photons allow an 
estimate of the muon energy. However, the resolution is limited by fluctuations of these 
processes. This estimate is a lower bound on the neutrino energy, because it is based on 
the muon energy at the detector. The interaction vertex may be far outside the detector. 

The v e and v T channels are different. The electron from a v e will generate an electro- 
magnetic cascade, which is confined to a volume of a few cubic meters. This cascade 
coincides with the hadronic cascade X of the primary interaction vertex. The optical 
signature is an expanding spherical shell of Cherenkov photons with a larger intensity 
in the forward direction. The tau from a v T will decay immediately and also generate 
a cascade. However, at energies > 1 PeV this cascade and the vertex are separated by 
several tens of meters, connected by a single track. This signature of two extremely bright 
cascades is unique for high-energy u T , and it is called a double bang event [3]. 

The measurement of cascade-like events is restricted to interactions close to the detector, 
thus requiring larger instrumented volumes than for detection. Also the accuracy of the 
direction measurement is worse for cascades than for long muon tracks. However, when 
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the flux is diffuse, the v e and v T channels also have clear advantages. The backgrounds 
from atmospheric neutrinos are smaller. The energy resolution is significantly better since 
the full energy is deposited in or near the detector. The cascade channel is sensitive to 
all neutrino flavors because the neutral current interactions also generate cascades. In 
this paper we focus on the reconstruction of muon tracks; details on the reconstruction of 
cascades are described in [4]. 

The most abundant events in AMANDA are atmospheric muons, created by cosmic rays 
interacting with the Earth's atmosphere. At the depth of AMANDA their rate exceeds the 
rate of muons from atmospheric neutrinos by five orders of magnitude. Since these muons 
are absorbed by the earth, a muon track from the lower hemisphere is a unique signature 
for a neutrino-induced muon 2 . The reconstruction procedure must have good angular 
resolution, good efficiency, and allow excellent rejection of down-going atmospheric muons. 

This paper describes the methods used to reconstruct muon tracks recorded in the AMAN- 
DA experiment. The AMANDA-II detector is introduced in section 2. The reconstruction 
algorithms and their implementation are described in sections 3 to 5. Section 6 summarizes 
event classes for which the reconstruction may fail and strategies to identify and eliminate 
such events. The performance of the reconstruction procedure is shown in section 7. We 
discuss possible improvements in section 8. 



2 The AMANDA Detector 

The AMANDA-II detector (see figure 2) has been operating since January 2000 with 677 
optical modules (OM) attached to 19 strings. Most of the OMs are located between 1500 m 
and 2000 m below the surface. Each OM is a glass pressure vessel, which contains an 8- 
inch hemispherical PMT and its electronics. AMANDA-B10 3 , the inner core of 302 OMs 
on 10 strings, has been operating since 1997. 

One unique feature of AMANDA is that it continuously measures atmospheric muons 
in coincidence with the South Pole Air Shower Experiment surface arrays SPASE-1 and 
SPASE-2 [7]. These muons are used to survey the detector and calibrate the angular 
resolution (see section 7 and [8,9]), while providing SPASE with additional information 
for cosmic ray composition studies [10]. 

The PMT signals are processed in a counting room at the surface of the ice. The analog 
signals are amplified and sent to a majority logic trigger [11]. There the pulses are discrim- 

2 Muon neutrinos above IPeV are absorbed by the Earth. At these ultra-high-energies (UHE), 
however, the the muon background from cosmic rays is small and UHE muons coming from 
the horizon and above are most likely created by UHE neutrinos. The search for these UHE 
neutrinos is described in [5,6]. 

3 Occasionally in the paper we will refer to this earlier detector instead of the full AMANDA-II 
detector. 
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Figure 2. The AMANDA-II detector. The scale is illustrated by the Eiffel tower at the left. 



inated and a trigger is formed if a minimum number of hit PMTs are observed within a 
time window of typically 2 /is. Typical trigger thresholds were 16 hit PMT for AMANDA- 
BIO and 24 for AMANDA-II. For each trigger the detector records the peak amplitude and 
up to 16 leading and trailing edge times for each discriminated signal. The time resolution 
achieved after calibration is at — 5 ns for the PMTs from the first 10 strings, which are 
read out via coaxial or twisted pair cables. For the remaining PMTs, which are read out 
with optical fibers the resolution is o t ~ 3.5 ns. In the cold environment of the deep ice 
the PMTs have low noise rates of typically 1 kHz. 

The timing and amplitude calibration, the array geometry, and the optical properties 
of the ice are determined by illuminating the array with known optical pulses from in 
situ sources [11]. Time offsets are also determined from the response to through-going 
atmospheric muons [12]. 

The optical absorption length in the ice is typically 110 m at 400 nm with a strong wave- 
length dependence. The effective scattering length at 400 nm is on average ~ 20 m. It 
is defined as A s /(1 — (cos# s )), where A s is the scattering length and 9 S is the scattering 
angle. The ice parameters vary strongly with depth due to horizontal ice layers, i.e., vari- 
ations in the concentration of impurities which reflect past geological events and climate 
changes [13-19]. 
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3 Reconstruction Algorithms 



The muon track reconstruction algorithm is a maximum likelihood procedure. Prior to 
reconstruction simple pattern recognition algorithms, discussed in section 4, generate the 
initial estimates required by the maximum likelihood reconstructions. 

3. 1 Likelihood Description 

The reconstruction of an event can be generalized to the problem of estimating a set of 
unknown parameters {a}, e.g. track parameters, given a set of experimentally measured 
values {x}. The parameters, {a}, are determined by maximizing the likelihood £(x|a) 
which for independent components Xi of x reduces to 

£(x|a) = IJpNa) , (2) 

i 

where p(xi\a) is the probability density function (p.d.f.) of observing the measured value 
xi for given values of the parameters {a} [20]. 




Figure 3. Cherenkov light front: definition of variables 

To simplify the discussion we assume that the Cherenkov radiation is generated by a 
single infinitely long muon track (with (3=1) and forms a cone. It is described by the 
following parameters: 

a = (r , t , p, E ) (3) 

and illustrated in figure 3. Here, r is an arbitrary point on the track. At time to, the 
muon passes r with energy E along a direction p. The geometrical coordinates contain 
five degrees of freedom. Along this track, Cherenkov photons are emitted at a fixed angle 
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6 C relative to p. Within the reconstruction algorithm it is possible to use a different 
coordinate system, e.g. a = (d, r], ...). The reconstruction is performed by minimizing 
— log(£) with respect to a. 

The values {x} presently recorded by AMANDA are the time ti and duration TOTi {Time 
Over Threshold) of each PMT signal, as well as the peak amplitude Ai of the largest 
pulse in each PMT. PMTs with no signal above threshold are also accounted for in the 
likelihood function. The hit times give the most relevant information. Therefore we will 
first concentrate on p(t|a). 



3.1.1 Time Likelihood 



According to the geometry in figure 3, photons are expected to arrive at OM % (at v\) at 
time 

p-(ri-r ) + dtan# c 
£ g co — to H , (4 J 

Cvac 

with c vac the vacuum speed of light 4 . It is convenient to define a relative arrival time, or 
time residual 

Ires = ^hit ^geo j (5) 

which is the difference between the observed hit time and the hit time expected for a 
"direct photon", a Cherenkov photon that travels undelayed directly from the muon to 
an OM without scattering. 
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Figure 4. Schematic distributions of arrival times £tes f° r different cases: Top left: PMT jitter. 
Top right: the effect of jitter and random noise. Bottom left: The effect of jitter and secondary 
cascades along the muon track. Bottom right: The effect of jitter and scattering. 



4 We note that equation 4 neglects the effect that Cherenkov light propagates with group 
velocity as pointed out in [21]. It was shown in [14] that for AMANDA this approximation is 
justified. 
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In the ideal case, the distribution, p(t res \a), would be a delta function. However, in a 
realistic experimental situation this distribution is broadened and distorted by several 
effects, which are illustrated in figure 4. The PMT jitter limits the timing resolution a t . 
Noise, e.g. dark noise of the PMT, leads to additional hits which are random in time. 
These effects can generate negative t res values, which would mimic unphysical causality 
violations. Secondary radiative energy losses along the muon trajectory create photons 
that arrive after the ideal Cherenkov cone. These processes are stochastic, and their 
relative photon yield fluctuates. 

In AMANDA, the dominant effect on photon arrival times is scattering in the ice 5 . The 
effect of scattering depends strongly on the distance, d, of the OM from the track as 
illustrated in figure 4. Since the PMTs have a non-uniform angular response, p(t re s) also 
depends on the orientation, 77, of the OM relative to the muon track (see figure 3). OMs 
facing away from the track can only see light that scatters back towards the PMT face. 
On average this effect shifts t rcs to later times and modifies the probability of a hit. 

The simplest time likelihood function is based on a likelihood constructed from p±, the 
p.d.f. for arrival times of single photons % at the locations of the hit OMs 

Whits 

Aime = II Pl(*res,i|a = di,T)i, . . .) . (6) 
i=l 

Note that one OM may contribute to this product with several hits. The function pi(t res ,i\ a) 
is obtained from the simulation of photon propagation through ice (see section 3.2). How- 
ever, this description is limited, because the electrical and optical signal channels can only 
resolve multiple photons separated by a few 100 ns and ~ 10 ns, respectively. Within this 
time window, only the arrival time of the first pulse is recorded. 

This first photon is usually less scattered than the average single photon, which modifies 
the probability distribution of the detected hit time. The arrival time distribution of the 
first of N photons is given by 

(roc \ (N— 1) 

/ Pl {t)dt) = A^ 1 (t tres )-(l-Pl(trcs)) (JV ~ 1) (7) 

Pi is the cumulative distribution of the single photon p.d.f.. The function pjv(t re s) is called 
the multi-photo-electron (MPE) p.d.f. and correspondingly defines £mpe- 

This concept can be extended to the more general case of p%(t Tes ), the p.d.f. for the k th 
photon out of a total of N to arrive at t res , given by 

P k N (tre S ) = N ■ ( N k Zl) * • " P ^)) {N ~ k) ' ( A^es))^ (8) 

PN^res) specifies the likelihood of arrival times of individual photoelectrons for averaged 
5 In water detectors this effect is neglected [22] or treated as a small correction [23]. 
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time series of N photoelectrons. With waveform recording the arrival times and amplitudes 
of individual pulses can be resolved. 

When the number of photoelectrons, N, is not measured precisely enough, multi-photon 
information can be included via another method. Instead of measuring N, the p.d.f. of 
the first photoelectron can be calculated by convolving the MPE p.d.f. Px(t,d) with the 
Poisson probability P^ olsson (//), where \x is the mean expected number of photoelectrons 
as a function of the distance, d. 

Pl(*») = = T^-PiM-e-M*™) (9) 

This result is called the Poisson Saturated Amplitude (PSA) p.d.f. [24, 25] and corre- 
spondingly defines £psa- The constant N — 1 — renormalizes the p.d.f. to unity. 

The probability of (uncorrelated) noise hits is small. They are further suppressed by a hit 
cleaning procedure (section 5.3), which is applied before reconstruction. They are included 
in the likelihood function by simply adding a constant p.d.f. po. 



3.1.2 Hit and No-Hit Likelihood 

The likelihood in the previous section relies only on the measured arrival times of photons. 
However, the topology of the hits is also important. PMTs with no hits near a hypothetical 
track or PMTs with hits far from the track are unlikely. 

A likelihood utilizing this information can be constructed as 

Cut = n phit,i • n p n °- hit '\ (io) 

i=l i=N ch +l 

where iV c h is the number of hit OMs and Aqm the number of operational OMs. The 
probabilities P hlt and p no ^ hlt of observing or not observing a hit depend on the track 
parameters a. Additional hits due to random noise are easily incorporated: p no_hlt — > 

pno— hit — pno— hit . pno— noise gj-^ phit ^ phit _ -y pno— hit 

Assuming that the probability is known for a single photon, the hit and no-hit 
probabilities of OMs for n photons can be calculated: 



and 



pno-Mt = (i - p^y 

(11) 

phit _ i pno— hit 

n n 

= l- (i-p^y . 

The number of photons, n, depends on E^, the energy of the muon: n = n(E^). For a 
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fixed track geometry, the likelihood (equation 10) can be used to reconstruct the muon 
energy. 



3.1.3 Amplitude Likelihood 

The peak amplitudes recorded by AMANDA can be fully incorporated in the likelihood [26], 
which is particularly useful for energy reconstruction. The likelihood can be written as 

W N ° M 

II W i ■ P i( A i) , (12) 

where Pi{AA is the probability that OM i observes an amplitude Ai, with A* = for 
unhit OMs. W and Wi are weight factors, which describe deviations of the individual OM 
and the total number of hit OMs from the expectation. Pj depends on the mean number 
/i of expected photoelectrons: 

P iiA) = P^. (l -P h .I^L. ( 13 ) 

The probability Pi{AA is normalized to the probability of observing the most likely ampli- 
tude (Ai). Pj th (/u) is the probability that a signal of \i does not produce a pulse amplitude 
above the discriminator threshold. As before, P hlt = 1 — p n °- hlt ; where the no-hit prob- 
ability is given by Poisson statistics: P no_hlt = exp(— fx) ■ (1 — p noisc ). The probability of 
Ai = is a special case: P(0) = P no_hlt _|_p hlt . p. th . Energy reconstructions based on this 
formulation of the likelihood will be referred to as Full E reco . 

An alternative energy reconstruction technique (see section 3.2.4) uses a neural net which 
is fed with energy sensitive parameters. 



3.1.4 Zenith Weighted (Bayesian) Likelihood 

Another extension of the likelihood [27-29] incorporates external information about the 
muon flux via Bayes' Theorem. This theorem states that for two hypotheses A and B, 

P{m= wm. (14) 

Identifying A with the track parameters a and B with the observations x, equation 14 gives 
the probability that the inferred muon track a was in fact responsible for the observed 
event x. P(x|a) is the probability that a, assumed to be true, would generate the event 
x — in other words, the likelihood described in the previous sections. -P(a) is the prior 
probability of observing the track a; i.e., the relative frequencies of different muon tracks 
as a function of their parameters. -P(x), which is independent of the track parameters a, 
is a normalization constant which ensures that equation 14 defines a proper probability. 
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Because the likelihood is only defined up to an arbitrary constant factor, this normalization 
may be ignored in the present context. 

In order to obtain P(a|x), one thus has to determine the prior probability distribution, 
-P(a), of how likely the various possible track directions are a priori. The reconstruction 
maximizes the product of the p.d.f. and the prior. 

The flux of muons deep underground is reasonably well known from previous experi- 
ments. Any point source of muons would be at most a small perturbation on the flux of 
penetrating atmospheric muons and muons created by atmospheric neutrinos. The most 
striking feature of the background flux from atmospheric muons is the strong dependence 
on zenith angle. For vertically down-going tracks it exceeds the flux from neutrino induced 
muons by about 5 orders of magnitude but becomes negligible for up-going tracks. This 
dependence, which is modeled by a Monte Carlo calculation [30], acts as a zenith depen- 
dent weight to the different muon hypotheses, a. With this particular choice, some tracks, 
which would otherwise reconstruct as up-going, reconstruct as down-going tracks. This 
greatly reduces the rate at which penetrating atmospheric muons are mis-reconstructed as 
up-going neutrino events [31]. In principle, a more accurate prior could be used. It would 
need to include the depth and energy dependence of the atmospheric muons as well as 
the angular dependence of atmospheric neutrino induced muons. 

Upon completion of this work, we learned that this technique was developed independently 
by the NEVOD neutrino detector collaboration [32] who were able to extract an atmo- 
spheric neutrino from a background of 10 10 atmospheric muons in a small (6 x 6 x 7.5m 3 ) 
surface detector. 



3.1.5 Combined Likelihoods 

The likelihood function Aimc of the hit times is the most important for track reconstruc- 
tion. However, it is useful to include other information like the hit probabilities. The 
combined p.d.f. from Equations 7 and 10 is 

^MPEeP Mt P no_Mt = ^MPE ' (Aiit) , (15) 

which is particularly effective. Here w is an optional weight factor which allows the ad- 
justment of the relative weight of the two likelihoods. This likelihood is sensitive not only 
to the track geometry but also to the energy of the muon. 

As discussed in section 3.1.4, the zenith angle dependent prior function, P(0), can be 
included as a multiplicative factor. This combination 

^Baycs = P(Q) ' Aime (16) 

has been used in the analysis of atmospheric neutrinos [30] . However, all of these improved 
likelihoods are limited by the underlying model assumption of a single muon track. 
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3.2 Likelihood Implementation 



The actual implementation of the likelihoods requires detailed knowledge of the photon 
propagation in the ice. On the other hand, efficiency considerations and numeric problems 
favor a simple and robust method. 

The photon hit probabilities and arrival time distributions are simulated as functions of 
all relevant parameters with a dedicated Monte Carlo simulation and archived in large 
look-up tables. This simulation is described in [26,30,33]. 

The AMANDA Collaboration has followed different strategies for incorporating this data 
into the reconstruction. In principle the probability density functions are taken directly 
from these archives. However, one has to face several technical difficulties due to the 
memory requirements of the archived tables, as well as numeric problems related to the 
normalization of interpolated bins and the calculation of multi-photon likelihoods. 

Alternatively, one can simplify the model and parametrize these archives with analyti- 
cal functions, which depend only on a reduced set of parameters. Comparisons of two 
independent parametrizations [24, 34] show that the direct and parametrized approaches 
yield similar results in terms of efficiency. This indicates that the parametrization itself is 
not limiting the reconstruction quality; rather, as mentioned earlier, the reconstruction is 
limited by the assumptions of the model being fit. Therefore, we will concentrate on only 
one parametrization. 



3.2.1 Analytical Parametrization 

A simple parametrization of the arrival time distributions can be achieved with the fol- 
lowing function, which we call Pandel function. It is a gamma distribution and its usage 
is motivated by an analysis of laser light signals in the BAIKAL experiment [35]. There, 
it was found that for the case of an isotropic, monochromatic and point-like light source, 
Pi (ires) can be expressed in the form 



(. / 1 , C me dium \ , d \ 
trcs • - + \ + T- 

Vr A a / X a (17) 

- N(d) T(d/X) [U) 

N(d) = e- d ' x * ■ (l + r " C " lcdium V" /A , (18) 
V A a / 

without special assumptions on the actual optical parameters. Here, c mcdium = c vac /n is 
the speed of light in ice, A a the absorption length, T(d/X) the Gamma function and N(d) a 
normalization factor, which is given by equation 18. This formulation has free parameters 
A and r, which are unspecified functions of the distance d and the other geometrical 
parameters. They are empirically determined by a Monte Carlo model. 
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The Pandel function has some convenient mathematical properties: it is normalized, it is 
easy to compute, and it can be integrated analytically over the time, t ies , which simplifies 
the construction of the multi-photon (MPE) time p.d.f.. For small distances the function 
has a pole at t = corresponding to a high probability of an unscattered photon. Going 
to larger values of d, longer delay times become more likely. For distances larger than the 
critical value d = A, the power index to t res changes sign, reflecting that the probability 
of undelayed photons vanishes: essentially all photons are delayed due to scattering. 

The large freedom in the choice of the two parameter functions r (units of time) and 
A (units of length) and the overall reasonable behavior is the motivation to use this 
function to parametrize not only the time p.d.f. for point-like sources, but also for muon 
tracks [34]. The Pandel function is fit to the distributions of delay times for fixed distances 
d and angles 77 (between the PMT axis and the Cherenkov cone). These distributions are 
previously obtained from a detailed photon propagation Monte Carlo for the Cherenkov 
light from muons. The free fit parameters are r, A, A a and the effective distance d e g, which 
will be introduced next. 




200 400 500 1000 1500 



time delay / ns time delay / ns 

Figure 5. Comparison of the parametrized Pandel function (dashed curves) with the detailed 
simulation (black histograms) at two distances d from the muon track. 

When investigating the fit results as a function of d and angle 77 (see figure 3), we observe 
that already for a simple ansatz of constant r, A and A a the optical properties in AMANDA 
are described sufficiently well within typical distances. The dependence on 77 is described 
by an effective distance d e fr which replaces d in equation 17. This means that the time 
delay distributions for backward illumination of the PMT is found to be similar to a 
head-on illumination at a larger distance. The following parameters are obtained for a 
specific ice model, and are currently used in the reconstruction: 

t = 557 ns d e ff — ao + a\ ■ d 

A = 33.3 m ai = 0.84 (19) 
A a = 98 m ao = 3.1 m — 3.9 m ■ cos(?7) + 4.6 m ■ cos 2 (t?) 

A comparison of the results from this parametrization with the full simulation is shown 
in figure 5 for two extreme distances. The simple approximation describes the behavior of 
the full simulation reasonably well. However, this simple overall description has a limited 
accuracy, especially for d ~ A (not shown). Reconstructions, based on the Pandel function 
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with different ice models, and a generic reconstruction, that uses the full simulation results, 
yield similar results. These comparisons indicate that the results of the reconstruction do 
not critically depend on the fine tuning of the underlying ice models, and it justifies the 
use of the above simple model. 



3.2.2 Extension to realistic PMT signals 

Although the Pandel function is the basis of a simple normalized likelihood, it has several 
deficiencies. It is not defined for negative t res , it ignores PMT jitter, and it has a pole at 
t res = 0, which causes numerical difficulties. These problems can be resolved by convolving 
the Pandel function, equation 17, with a Gaussian, which accounts for the PMT jitter. 
Unfortunately such convolution requires significant computing time. 

Instead the Pandel function is modified by extending it to negative times, t res < 0, with a 
(half) Gaussian of width o~ g . The effects of PMT jitter are only relevant for small values of 
t res . For times t TCS > ti the original function is used, and the two parts are connected by a 
spline interpolation (3rd order polynomial). The result, P(t TCS ), is called upandel function. 

Using ti = and requiring further a smooth interpolation and the normalization to 

be unchanged, the polynomial coefficients aj and normalization of the Gaussian N g can be 
calculated analytically [34]. The free parameter o~ g includes all timing uncertainties, not 
just the PMT jitter. Good reconstruction results are achieved for a large range 10 ns < 
o-g < 20 ns. 



3.2.3 p hlt p no - lut Parametnzation 



The normalization N(d) in equation 18 is used to construct a hit probability function, 
Phit- The function P^ 1 * with Pj ut = N(d), is fit to the hit probability determined by the full 
AMANDA detector simulation, as a function of distance, orientation and muon energy. The 
free parameters are the Pandel parameters r and A, A a , d and h. The effective distance d, 
is similar to the effective distance in the Pandel parametrization. We define n as the power 
index of equation 11, which corresponds to an effective number of photons. It is important 
to understand that N(d) is not a hit probability and h is not just a number of photons. 
They are constructs, that are calibrated with a Monte Carlo simulation. Technically, the 
power index, n = N in equation 11, factorizes into n(r],E^) = e^P^) • n(P^). The 
variable n, where n = n(E), is related to the number of photons incident on the PMT 
and its absolute efficiency The factor e(r], E^) is related to the orientation dependent 
PMT sensitivity but also accounts for the energy dependent angular emission profile of 
photons with respect to the bare muon. 
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3.2.4 Energy Reconstruction 

The reconstruction of the track geometry is a search for five parameters. If the muon 
energy is added as a fit parameter, the minimization is significantly slower. Therefore, 
the energy reconstruction is performed in two steps. First, the track geometry is recon- 
structed without the energy parameter. Then these geometric parameters are used in an 
energy reconstruction, that only determines the energy. However, if the time likelihood 
utilizes amplitude information, e.g. in the combined likelihood, equation 7, or the PSA 
likelihood, equation 9, the track parameters also depend on the energy. In this case the en- 
ergy and geometric parameters must be reconstructed together. Currently, three different 
approaches are used to reconstruct the muon energy. They are compared in section 7.3. 

(1) The simplest method utilizes the p hlt p n °- hlt reconstruction (see sections 3.1.2 and 
3.2.3). 

(2) The Full E reco method (see section 3.1.3) models the measured amplitudes in a like- 
lihood for reconstructing the energy. This algorithm performs better, but it is more 
dependent on the quality of the amplitude calibration of the OMs. 

(3) An alternative way to measure the energy is based on a neural network [36]. The 
neural network uses 6-6-3-1 and 6-3-5-1 feed-forward architecture for AMANDA-BIO 
and AMANDA-II, respectively. The energy correlated variables which are used as 
input are the mean of the measured amplitudes (ADC), the mean and RMS of the 
arrival times (LE) or pulse durations (TOT), the total number of signals, the number 
of OMs hit and the number of OMs with exactly one hit. 

Less challenging than a full reconstruction, a lower energy threshold is determined by 
requiring a minimum number of hit OMs. The number of hit OMs is correlated with 
the energy of the muon. Since celestial neutrinos are believed to have a substantially 
harder spectrum than atmospheric neutrinos, an excess of high multiplicity events would 
indicate that a hard celestial source exists. Values for this parameter determined from 
AMANDA data already set a tight upper limit on the diffuse flux of high-energy celestial 
neutrinos [37]. 



3.2.5 Cascade Reconstruction 

The reconstruction of cascade like events is described in detail elsewhere [4]. The basic 
approach is similar to the track reconstruction. It assumes events form a point light 
source with photons propagating spherically outside with a higher intensity in the forward 
direction. The cascade reconstruction also uses the Pandel function (see equation 17) with 
parameters that are specific for cascades. 

In several muon analyses, a cascade fit is used as a competing model. In cases where the 
cascade fit achieves a better likelihood than the track reconstruction, the track hypothesis 
is rejected. In particular this is used as a selection criterion to reject background events 
which are mis-reconstructed due to bright secondary cascades. 
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4 First Guess Pattern Recognition 



The likelihood reconstructions need an initial track hypothesis to start the minimization. 
The initial track is derived from first guess methods, which are fast analytic algorithms 
that do not require an initial track. 



4-1 Direct Walk 



A very efficient first guess method is the direct walk algorithm. It is a pattern recognition 
algorithm based on carefully selected hits, which were most likely caused by direct photons. 

The four step procedure starts by selecting track elements, the straight line between any 
two hit OMs at distance d, which are hit with a time difference 

|At| < — + 30ns with d>50m. (20) 

Cyac 

The known positions of the OMs define the track element direction (0,<f>). The vertex 
position (x, y, z) is taken at the center between the two OMs. The time at the vertex t 
is defined as the average of the two hit times. 

In a next step, the number of associated hits (AH) are calculated for each track element. 
Associated hits are those with —30 ns < t rC s < 300 ns and d < 25 m- (t rcs + 30) 1 / 4 (t in ns), 
where d is the distance between hit OM and track element and t res is the time residual, 
which is defined in Equation 5. After selecting these associated hits, track elements of poor 

quality are rejected by requiring: Nab > 10 and a L = ^( ]\?ah ~ (L)) 2 ) > 20 m. 

Here, the "lever arm" Lj is the distance between the vertex of the track element and the 
point on the track element which is closest to OM % and (L) is the average of all Lj-values. 
Track elements that fulfill these criteria qualify as track candidates (TC). 

Frequently, more than one track candidate is found. In this cluster search is 

performed for all track candidates that fulfill the quality criterion: 



QTC - °- 7 ' Qmax , where 

Qmax = max(Q^Q) and (21) 
= min(Npjj[, 0.3 m" 1 • a L + 7). 

In the cluster search, the "neighbors" of each track candidate are counted, where neighbors 
are track candidates with space angle differences of less than 15°. The cluster with the 
largest number of track candidates is selected. 
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In the final step, the average direction of all track candidates inside the cluster defines 
the initial track direction. The track vertex and time are taken from the central track 
candidate in the cluster. Well separated clusters can be used to identify independent 
muon tracks in events which contain multiple muons (see section 6.1). 



4.2 Line-Fit 



The line-fit [38] algorithm produces an initial track on the basis of the hit times with 
an optional amplitude weight. It ignores the geometry of the Cherenkov cone and the 
optical properties of the medium and assumes light traveling with a velocity v along a 
1-dimensional path through the detector. The locations of each PMT, r*, which are hit at 
a time tj can be connected by a line: 

r,^ r + v • U , (22) 

A x 2 to be minimized is defined as: 

iv hit 

^^(r.-r-vt,) 2 , (23) 
i=i 

where A^it is the number of hits. The x 2 is minimized by differentiation with respect to 
the free fit parameters r and v. This can be solved analytically: 

r = <r<> - v ■ ( ti ) and v= ^'^^fy^ , (24) 

where (xi) = Xi denotes the mean of parameter x with respect to all hits. 

The line-fit thus yields a vertex point r, and a direction e = Vlf/|vlf|- The zenith angle 
is given by #lf = — arccos(t> 2 /|v L F|)- 

The time residuals (equation 5) for this initial track generally do not follow the distribution 
expected for a Cherenkov model. If the to parameter of the initial track is shifted to 
better agree with a Cherenkov model, subsequent reconstructions converge better (see 
section 5.2.3). 

The absolute speed i>lf = |v|, of the line-fit is the mean speed of the light propagating 
through the 1-dimensional detector projection. Spherical events (cascades) and high en- 
ergy muons have low i>lf values, and thin, long events (minimally ionizing muon tracks) 
have large values. 
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4-3 Dipole Algorithm 



The dipole algorithm considers the unit vector from one hit OM to the subsequently hit 
OM as an individual dipole moment. Averaging over all individual dipole moments yields 
the global moment M. It is calculated in two steps. First, all hits are sorted according to 
their hit times. Then a dipole-moment M is calculated: 

M = — • V * i ~\ . (25) 
N ch -1 felri-ri-il v ; 

It can be expressed via an absolute value Mda = ]M| and two angles #da and 0da- These 
angles define the initial track. 

The dipole algorithm does not generate as good an initial track as the direct walk or the 
line-fit, but it is less vulnerable to a specific class of background events: almost coincident 
atmospheric muons from independent air showers in which the first muon hits the bottom 
and the second muon hits the top of the detector. 



4-4 Inertia Tensor Algorithm 



The inertia tensor algorithm is based on a mechanical picture. The pulse amplitude from 
a PMT at corresponds to a virtual mass dj at r^. One can then define the tensor of 
inertia I of that virtual mass distribution. The origin is the center of gravity (COG) of 
the mass distribution. The COG-coordinates and the tensor of inertia components are 
given by: 

COG ee £ ( ai ) w ■ Vi and 

Zl ( 26 ) 
= E(ai) w ■[5 kl -(r t ) 2 -rf -4 . 
i=i 

The amplitude weight w > can be chosen arbitrarily. The most common settings are 
w = (ignoring the amplitudes) and w = 1 (setting the virtual masses equal to the 
amplitudes). The tensor of inertia has three eigenvalues Ij, je{l,2,3}, corresponding to 
its three main axes ej. The smallest eigenvalue I\ corresponds to the longest axis ex. In 
the case of a long track-like event I\ <^ {I 2 , h} and ei approximates the direction of 
the track. The ambiguity in the direction along the e\ axis is resolved by choosing the 
direction where the average OM hit time is latest. In the case of a cascade-like event, 
I ± « J 2 pa I 3 . The ratios between the Ij can be used to determine the sphericity of the 
event. 
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5 Aspects of the technical Implementation 



5.1 Reconstruction Framework 
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Figure 6. Schematic principle of the reconstruction chain 

The basic reconstruction procedure, sketched in figure 6, is sequential. A fast reconstruc- 
tion program calculates the initial track hypothesis for the likelihood reconstruction. All 
reconstruction programs may use a reduced set of hits in order to suppress noise hits and 
other detector artifacts. Event selection criteria can be applied after each step to reduce 
the event sample, and allow more time consuming calculations at later reconstruction 
stages. This procedure may iterate with more sophisticated but slower algorithms ana- 
lyzing previous results. The final step is usually the production of Data Summary Tape 
(DST) like information, usually in form of PAW N-tuples [39]. A detailed description of 
this procedure can be found in [40]. 

The reconstruction framework is implemented with the recoos program [41], which is based 
on the rdmc library [42] and the SiEGMuND software package [43]. The recoos program 
is highly modular, which allows flexibility in the choice and combination of algorithms. 



5.2 Likelihood Maximization 



The aim of the reconstruction is to find the track hypothesis which corresponds to the 
maximum likelihood. This is done by minimizing — log(£) with respect to the track pa- 
rameters. We have implemented several minimization procedures. 

The likelihood space for AMANDA events is often characterized by several minima. Local 
likelihood minima can arise due to symmetries in the detector, especially in the azimuth 
angle, or due to unexpected hit times caused by scattering. In the example, shown in 
figure 7, the reconstruction converged on a local minimum, because of non-optimal starting 
values. Several techniques, which are used to find the global minimum, are here presented. 
In particular, the iterative reconstruction, section 5.2.2, solves the problem and converges 
to the global minimum. One generally assumes that the global minimum corresponds to 
the true solution, but this is not always correct due to stochastic nature of light emission 
and detection. Such events cannot be reconstructed properly and have to be rejected using 
quality parameters (see section 6.2). 
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Figure 7. An example of the likelihood space (1-dimensional projection) for a specific AMANDA 
event. Shown is — log(£) as function of the zenith angle. Each point represents a fit, for which 
the zenith angle was fixed and the other track parameters were allowed to vary in order to find 
the best minimum. A local minimum which was found by a gradient likelihood minimization is 
indicated by a fitted parabola. Improved methods that avoid this are described in the text. 

5.2.1 Minimization algorithms 

The reconstruction framework allows us to use and compare these numerical minimization 
algorithms: Simplex [44], Powell's [44], Minuit [45] (using the minimize method), and 
Simulated annealing [44]. The Simplex algorithm is the fastest algorithm. Powell's method 
and Minuit are ~5 times slower than the Simplex algorithm. The reconstruction results 
from Minuit and the Simplex algorithm are nearly identical and almost as good as the 
Powell results. Exceptions occur in less than 1% of the cases, when these methods fail 
and stop at the extreme zenith angles 6 = 0° and 9 = 180°. The Simulated annealing 
algorithm is less sensitive to local minima than the other algorithms, but it is much 
slower and requires fine-tuning. 



5.2.2 Iterative Reconstruction 

The iterative reconstruction algorithm successfully copes with the problem of local minima 
and extreme zenith angles by performing multiple reconstructions of the same event. Each 
reconstruction starts with a different initial track. Therefore, the fast Simplex algorithm 
is sufficient. 

The ability to find the global minimum depends strongly on the quality of the initial track. 
A systematic scan of the full parameter space for initial seeds is not feasible. Instead the 
iterative algorithm concentrates on the direction angles, zenith and azimuth, and uses 
reasonable values for the spatial coordinates. The following procedure yields good results. 

The result of a first minimization is saved as a reference. Then both direction angles are 
randomly selected. The track point, r , is transformed to the point on the new track, 
which is closest to the center of gravity of hits. The time, to, of this point is shifted to 
match the Cherenkov expectation (see section 5.2.3). Then a new minimization is started. 
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If the minimum is less than the reference minimum, it is saved as the new reference. This 
procedure is iterated n times, and the best minima found for zenith angles above and 
below the horizon, are saved, and used to generate an important selection parameter (see 
section 6.2). 

This algorithm substantially reduces the number of false minima found, after a few itera- 
tions. For n = 6 roughly 95% of the results are in the vicinity of the asymptotic optimum 
for n — > oo. For n ~ 20 more than 99% of the results are the global minimum. Despite the 
fast convergence, the iterative reconstruction requires significant CPU time, which limits 
its use to reduced data sets. 

5.2.3 Lateral shift and time residual 

The efficiency of finding the global minimum of the likelihood function can be improved 
by translating the arbitrary vertex and/or time origin of the output track from the first 
guess algorithm before application of the full maximum likelihood method. 

Transformation of r . In general this "vertex" point is arbitrary in the infinite track 
approximation used, and first guess methods may produce positions distant from the 
detector. During the likelihood minimization, numerical errors can be avoided and the 
convergence improved by shifting this point along the direction of the track towards the 
point closest to the center of gravity of hits (see equation 26). The vertex time to is 
transformed accordingly: A(t ) = A(r )/c. 

Transformation of to- The time to obtained from first guess algorithms is not calculated 
from a full Cherenkov model. The efficiency of the likelihood reconstruction can be im- 
proved by shifting the t such that the time residuals, equation 5, fit better to a Cherenkov 
model. In particular it is useful to avoid negative t TCS , which would correspond to causal- 
ity violations. This can be achieved by transforming t — > t — £~ s , where t~ s is the most 
negative time residual. 

5.2.4 Coordinates and restricted parameters 

The track coordinates a, which are used by the likelihood, are independent of the coor- 
dinates actually chosen for the minimization. Therefore, the coordinate system can be 
chosen arbitrarily. Any of the parameters in this coordinate system can be kept fixed. 
During the minimization parameterization functions translate the coordinates as neces- 
sary. The most commonly used coordinates are ro and the zenith and azimuth angles 9, 
<P. 

The freedom in the choice of coordinates can be used to improve the numerical minimiza- 
tion, for systematic studies, or to fix certain parameters according to external knowledge. 
An example is the reconstruction of coincident events with the SPASE surface arrays [8]. 
Here, we fix the location of the trajectory to coincide with the core location at the surface 
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as measured by SPASE. Then, the direction is determined with the AMANDA reconstruc- 
tion subject to this constraint. 

Under certain circumstances the allowed range of the reconstruction parameters is re- 
stricted. The most important example here is to restrict the reconstructed zenith angle 
to above or below the horizon, to find the most likely up- or down-going tracks, respec- 
tively. Comparing the quality of the two solutions can be used for background rejection. 
Technically the constrained fit is accomplished by multiplying the likelihood by a prior, 
which is zero outside the allowed parameter range. 



5.3 Preprocessing and Hit cleaning 

The data must be filtered and calibrated before reconstruction. Defective OMs are re- 
moved, and the amplitudes and hit times are calibrated. A hit cleaning procedure iden- 
tifies and flags hits which appear to be noise or electronic effects, such as cross talk or 
after-pulsing. These hits are not used in the reconstruction, but they are retained for 
post-reconstruction analysis. 

The hit cleaning procedure can be based on simple and robust algorithms, because the 
PMTs have low noise rates. Noise and after-pulse hits are strongly suppressed by rejecting 
hits that are isolated in time and space from other signals in the detector. Typically a hit 
is considered to be noise if there is no hit within a distance of 60 m to 100 m and a time of 
±300 ns to ±600 ns. Cross talk hits are identified by examining the amplitudes and pulse 
widths of the individual pulses and by analyzing the correlations of uncalibrated hit times 
with hits of large amplitude in channel combinations which are known to cross talk to 
each other. The required cross talk correlation map was determined independently in a 
dedicated calibration campaign. 



5.4 Processing Speeds 

The first guess algorithms are sufficiently fast that the execution time is dominated by file 
input/output and the software framework. Typical fit times are ~ 20 ms per event on a 
850 MHz Pentium-Ill Linux PC. The processing speed of the likelihood reconstructions can 
vary significantly depending on the number of free parameters, the number of iterations, 
the minimization algorithm, and the experimental parameters like the number of hit OMs. 
These effects dominate the differences in processing speeds due to different reconstruction 
algorithms. The typical execution time for a 16-fold iterative likelihood reconstructions 
using the simplex minimizer to reconstruct the 5 free track parameters is ~ 250 ms per 
event. 
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6 Background Rejection 



The performance of the reconstruction depends strongly on the quality and background 
selection criteria. The major classes of background events in AMANDA (see section 6.1) are 
suppressed by the quality parameters presented in section 6.2. Optimization strategies for 
the selection criteria are summarized in section 6.3. Finally, we evaluate the reconstruction 
performance in section 7. 



6.1 Background Classes 

Most background events from atmospheric muons are well reconstructed and can be re- 
jected by selecting up-going reconstructed events. However, there is a small fraction of 
mis-reconstructed events, amounting to about 1CT 2 for the unbiased and about 10" 4 for 
the zenith-weighted reconstruction. These events are rejected by additional selection cri- 
teria described in section 6.2. These background events are classified as follows. 

Nearly horizontal muons: These events have true incident angles close to the horizon. 
A small error in the reconstruction causes them to appear as up-going. These events 
are not severely mis-reconstructed, but occur due to the finite angular resolution. 

Muon bundles: The spatial separation between multiple muons from a single air shower, 
a muon bundle, is usually small enough that the event can be described by a single bright 
muon track. If the separation is too large, the reconstruction fails. 

Cascades: Bright stochastic energy losses (e.g. bremsstrahlung) produce additional light, 
which distorts the Cherenkov cone from the muon. Cascades emit most of their light 
with the same Cherenkov angle as the muon, but some light is emitted at other an- 
gles. These secondary events can cause the reconstruction to fail, especially when the 
cascade(s) produce more light than the muon itself. A special class of these events are 
muons which pass outside the detector and release a bright cascade, which can mimic 
an up-going hit pattern. 

Stopping muons: Over the depth of the detector the muon flux changes by a factor of 
~2, since muons lose their energy and stop. These muons can create an up-going hit 
pattern, especially when the muon stops just before entering the detector from the side. 

Scattering layers: The scattering of light in the polar ice cap varies with depth. Light 
from bright events, can mimic an up-going hit pattern, in particular when it traverses 
layers of higher scattering. 

Corner clippers: These are events where the muon passes diagonally below the detector. 
The light travel upwards through the detector mimicking an up-going muon. 

Uncorrelated coincident muons: Due to the large size of the AMANDA detector, the 
probability of muons from two independent air showers forming a single event is small 
on the trigger level but not negligible. If an initial muon traverses the bottom of the 
detector and a later muon traverses the top, the combination can be reconstructed as 
an up-going muon. 
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Electronic artifacts: Noise, cross talk and other transient electronic malfunctions are 
generally small effects, but they can occasionally produce hits, which distort the time 
pattern. Such effects become important after a selection process of several orders of 
magnitude. 



6.2 Quality Parameters 



Background events, which pass a zenith angle selection, need to be rejected by applying 
selection criteria on quality parameters. These parameters usually evaluate information, 
which is not optimally exploited in the reconstruction. The detailed choice of quality 
parameters is specific to each analysis. Here, we summarize the most important categories. 

The number of direct hits, N^ti : £2), is the number of hits with small time residuals: 
t\ < t ms < t 2 (see equation 5). Un-scattered photons provide the best information for the 
reconstruction, and a large number of indicates high quality information in the event. 
Empirically reasonable values are t\ ~ —15 ns and ti between = +25 ns and +150 ns, 
depending on the specific analysis. 

The length of the event L is obtained by projecting each hit OM onto the reconstructed 
track and taking the distance between the two outermost of these points. L can be con- 
sidered as the "lever arm" of the reconstruction. Larger values corresponding to a more 
robust and precise reconstruction of the track's direction. This parameter is particularly 
powerful when calculated for direct hits only, and is then referred to as Ldir(ii '■ h)- Length 
requirements are efficient against corner clippers, stopping muons and cascades. 

The absolute value of the likelihood at the maximum is a good parameter to evaluate the 
quality of a reconstruction. Here, a useful observable is the likelihood parameter L which 
is defined as 

LS -^, (27) 

1 ' tree 

where iVf ree is the degrees of freedom (e.g. 7Vf ree = A^its — 5 for a track reconstruction). For 
Gaussian probability distributions this expression corresponds to the reduced chi-square. 
L can be used as a selection parameter, smaller values corresponding to higher quality. A 
selection of events with good |_^ hlt -P no hlt values is efficient against stopping muons. 

Comparing L from different reconstructions is a powerful technique. Cascade-like events 
will have a better likelihood from a cascade reconstruction than one from a track recon- 
struction. 

Another efficient rejection method is to compare L for the best up-going versus the best 
down-going reconstruction of a single event. If the up-going reconstruction is not signifi- 
cantly better than the down-going reconstruction, the event is rejected. These values can 
be obtained from the iterative reconstruction method (section 5.2.2) or by restricting the 
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parameter space. This method is particularly powerful when the down-going reconstruc- 
tion uses a zenith weighted likelihood (section 3.1.4). 



The reconstruction methods consider the p.d.f. for each hit separately but ignore corre- 
lations. Therefore, the reconstructions assign the same likelihood to tracks where all hits 
cluster at one end of the reconstructed track and tracks where the same number of hits 
are smoothly distributed along the track. The latter hit pattern indicates a successful 
track reconstruction, while the former hit pattern may be caused by a background event. 
The smoothness parameter S was inspired by the Kolmogorov-Smirnov test of the consis- 
tency of two distributions. S is a measure of the consistency of the observed hit pattern 
with the hypothesis of constant light emission by a muon. The simplest definition of the 
smoothness S is S — SJ 1 ^, where S™ ax is that Sj, which has the largest absolute value, 
and Sj is defined as 

Si = AjZ ^~ ~ r- ■ (28) 
3 N-l l N K J 

lj is the distance along the track between the points of closest approach of the track to 

the first and the j th hit module, with the hits taken in order of their projected position 

on the track. N is the total number of hits. Tracks with hits clustered at the beginning or 

end of the track have S approaching +1 or —1, respectively. High quality tracks with S 

close to zero, have hits equally spaced along the track. A graphical representation of the 

smoothness construction can be found in [30]. 

Extensions of this smoothness parameter include the restriction of the calculation to direct 
hits only or using the distribution of hit times U instead of the distances /«. 

A particularly important extension is S p * . In order to account for the granularity and 
asymmetric geometry of the detector one can replace the above formulation with one 
that models the hit smoothness expectation for the actual geometry of the assumed muon 
track. This can be accomplished by using the hit probabilities of all Nqm, the number of 
operational OMs, (ordered along the track) as weights: S phlt = max(S' phIt ) with 



2 ■ ° M A- ^ °j M _phit,i v ' 

At = 1, if the OM % was hit and otherwise and Phit,i is the probability for OM % to 
be hit given the reconstructed track. The hit probabilities are calculated according to 
the algorithm in section 3.2.3. Smoothness selections are very efficient against secondary 
cascades, stopping muons and coincident muons from independent air showers. 

Interesting AMANDA events are analyzed with multiple reconstruction algorithms. An 
event is most likely to have been reconstructed correctly, if the different algorithms produce 
consistent results. 

For two reconstructions with directions ei and e 2 , the space angle between them is given by 
\I/ = arccos (ei • e 2 ), which should be reasonably small for successful reconstructions. This 
concept can be extended to multiple reconstructions and their angular deviations from the 
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average direction. For n different reconstructed directions, e^, the average reconstructed 
direction, E, is given by E = J2? e i/\ Y% G i\- We can define the parameter 



VPi describes the average space angle between the individual reconstructions and E. 
is a different parameter, which treats the deviations between E and the as "errors" 
and adds them quadratically. Small values of ^ or ^ indicate consistent reconstruction 
results. 

The \1/ parameters are a mathematically correct consistency check only when comparing 
the results of uncorrelated reconstructions of the same intrinsic resolutions. This is not 
the case when comparing different AMANDA reconstructions. Irrespective of the validity 
of such an interpretation, ^ or \l/ 2 are very efficient selection criteria, especially against 
almost horizontal muons and wide muon bundles. 

A few additional selection parameters are closely related to first guess methods. The ratio 
of the eigenvalues of the tensor of inertia (see section 4.4) are a measure of the sphericity of 
the event topology, which is an efficient selection parameter against cascade backgrounds. 
Tracks reconstructed as down-going by the dipole fit (see section 4.3) that have a non- 
negligible dipole moment, M^a = indicate coincident muons from independent air 
showers. Larger values of the line-fit speed vlf (see section 4.2) are an indication for longer 
muon-like, smaller values for more spherical cascade-like events. 

Finally, two approaches evaluate the "intrinsic resolution" or "stability" of the recon- 
struction of each event. One approach quantifies the sharpness of the minimum found by 
the minimizers in — log(£) by fitting a paraboloid to it. The fitted parameters can then 
be used to classify the sharpness of the minimum. The other approach splits an event 
into sub-events (for example, containing odd- vs. even-numbered hits) and reconstructs 
the sub-events. If the reconstructed directions of the sub-events are different, then the 
reconstruction of the full event has a larger uncertainty. 

6.3 Analysis Strategies 

Analyses that search for neutrino induced muons must cope with a large background 
of atmospheric muons. The optimal choice of reconstruction and selection criteria varies 
strongly with different expectations for the energy and angular distribution of the signal 
events. The goal is to optimize the signal efficiency over the background or noise (square 
root of the background) based on sets of signal and background data. 

• The selection criteria for background sensitive variables may be adjusted individually 
such that a specified fraction of signal events pass. After these first level criteria are 
set, the adjustment is repeated until the desired background rejection is reached. Each 




(30) 
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iteration defines a "cut-level" , which corresponds to data sets of increasing purity. This 
simple method is used to derive a defined set of selection parameters for the performance 
section 7. However, less efficient criteria are mixed with more efficient criteria, and 
correlations of the variables are not taken into account. Therefore, this method does 
not achieve the optimum signal efficiency. 

An improvement to this method has been demonstrated in an AMANDA point source 
analysis [46,47]. Here, a selection criterion is only applied to the most sensitive variable, 
and the most sensitive variable is determined at each cut level. An interesting aspect 
of this point source search is that the experimental data themselves can be used as 
a background sample, which reduces systematic uncertainties from the background 
simulation. The selection criteria are not optimized with respect to signal purity but 
with respect to an optimal significance of a possible signal. 

Another approach is to combine the selection parameters into a single selection param- 
eter, called event quality. This can be done by rescaling and normalizing each of the 
selection parameters according to the cumulative distribution of the signal expectation. 
The AMANDA analyses of atmospheric neutrinos [30,48] used this technique. 
Additional approaches use discriminant analysis [49] or neural nets [6,50,51] to opti- 
mize the efficiency while taking into account the correlations between selection criteria 
and their individual selectivity. However, both methods depend critically on a good 
agreement between experiment and simulation. These methods quantify the efficiency 
of each parameter by including and excluding it from the optimization procedure. 
The "CutEval" method finds the optimum combination of selection parameters and 
cut values by numerically maximizing a significance function, Q. An example is Q = 
S/ VB, where S is the number of signal events, and B is the number of background 
events after selection. 

The implementation proceeds in several steps. First, the most efficient selection pa- 
rameter, Ci, is the parameter that individually maximizes Q. The next parameter, C2, 
is the parameter that maximizes Q in conjunction with C\. More parameters are suc- 
cessively determined until the addition of a new parameter fails to improve Q. This 
procedure takes correlations between the selection criteria into account. The final num- 
ber of selection parameters is reduced to a minimum, while maximizing the efficiency. 
Next, the optimal selection for this combination parameters is computed as a function 
of a boundary condition (e.g. the maximum number of accepted background events). 
This boundary condition is also used to define a single quality parameter. 

Such a formalized procedure has to be carefully monitored, e.g. to handle potentially 
un-simulated experimental effects. The CutEval procedure is monitored by defining 
different, complementary optimization functions Q, which allow real and simulated data 
to be compared [30,52-55]. 
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7 Performance 



This section describes the performance of the reconstruction methods. It is based on 
illustrative data selections, and the actual performance of a dedicated analysis can be 
different. Unless noted otherwise, the data shown is from Monte Carlo simulations of 
atmospheric neutrinos for AMANDA-II. 

1.1 First Guess Algorithms 

Since the first guess algorithms are used as a starting point for the full reconstruction, 
they should provide a reasonable estimate of the track coordinates. Also, these algorithms 
are used as the basis of early level filtering, and therefore need to be sufficiently accurate 
for that purpose, i.e. they should at least reconstruct the events in the correct hemisphere. 



reconstruction 


atm. fj, 


atm. v 


direct walk 


1.5% 


93% 


line- fit 


4.8% 


85% 


dipole algorithm 


16.8% 


78% 



Table 1 

The atmospheric muon and atmospheric neutrino detection efficiencies for a selection at 9 > 80° 
for the first guess algorithms. 

As an example, Table 1 gives the passing efficiencies with respect to the AMANDA-II 
trigger for atmospheric neutrinos (signal) and atmospheric muons (background) for the 
first guess methods (see section 4), after the selection of events with calculated zenith 
angles larger 80°. The direct walk algorithm gives the best background suppression and 
the highest atmospheric neutrino passing rate. Correspondingly, it also gives the best 
initial tracks to the likelihood reconstructions. 

7.2 Pointing accuracy of the Track Reconstruction 

The angular accuracy of the reconstruction can be expressed in terms of a point spread 
function, which is given by the space angle deviation \1/ between the true and the re- 
constructed direction of a muon corrected for solid angle. The space angle deviation is a 
combined result of two effects: a systematic shift in the direction and a random spread 
around this shift. In a point source analysis, for example, it is possible to correct for 
systematic shifts and be limited by the point spread function alone [47]. 

The zenith and space angular deviations are shown in figures 8 and 9. They are obtained by 
the reconstruction algorithms as used in AMANDA-BIO. The same event selection is used 
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Figure 8. The zenith angle deviations for various reconstructions of AMANDA-BIO. The result 
of an atmospheric neutrino simulation after the selection criteria of [30] is shown. The fits are a 
line- fit (LF), an iterated upandel fit (LH), an iterated zenith- weighted upandel fit and a MPE 
fit. 
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Figure 9. The distribution of space angle deviations for various reconstructions of AMANDA-BIO. 
The result of an atmospheric neutrino simulation after the selection criteria of [30] is shown. 
The fits are a line- fit (LF), an iterated upandel fit (LH), an iterated zenith- weighted upandel fit 
and a MPE fit. 

for all. As a general observation, the distributions of deviations for different reconstruction 
algorithms is surprisingly similar after a particular selection. Larger differences are usually 
seen in the selection efficiencies. A similar behavior is observed for AMANDA-II. 



The dependence of the space angle deviation for the full AMANDA-II detector on the 
cut level 6 for the LH reconstruction is shown in figure 10. The tighter the selection 



The cut levels defined here are typical and intended as demonstrating example. We use typical 
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Figure 10. The dependence of the space angle deviation of the LH reconstruction in AMANDA-II 
on the event selection (cut levels). 

criteria, the better the angular resolution. The same general trend is true for the other 
reconstructions. Tight criteria select events with unambiguous hit topologies, which are 
reconstructed better. The results for cut level 6 are shown in figures 11-13 as function of 
the energy and the zenith angle. 



The angular resolution (see figure 11) has a weak energy dependence. The energy of the 
muon is taken at the point of its closest approach to the detector center. Best results are 
achieved for energies of 100 GeV to 10 TeV. At energies < 100 GeV, the muons have paths 
shorter than the full detector, which limits the angular resolution. At energies > 10 TeV, 
more light is emitted due to individual stochastic energy loss processes along the muon 
track. Here, the hit pattern is not correctly described by the underlying reconstruction 
assumption of a bare muon track (see section 6.1). 

The space angular resolution depends on the incident muon zenith angle (see figure 12). 
Again this is shown only for the LH reconstruction, the other reconstructions are similar. 
Up-going muons with cos 9^ ~ —0.7 are best reconstructed, and horizontal muons are 
the worst, because of the geometry of the AMANDA-II detector. Nearly vertical events 
with cos# M ~ —1 have a poorer angular resolution, because they illuminate fewer strings, 
which can cause ambiguities in the azimuth. 

Systematic shifts also degrade the angular resolution. AMANDA observes a small zenith 

selection parameters from section 6.2: the reconstructed zenith angle, # DW > 80°, # LH > 80°, 
iV ch , iV£f(-15 : 25), Ldir(~ 15 : 75 )> LLH > sLH and ^i(DW,LH,MPE). Our goal here is to 
illustrate the analysis, and we do not optimize with respect to efficiency and angular resolution. 
Instead each individual criterion is enforced in such a way that 95% of the events from the 
previous level would pass, and correlations between the parameters are ignored. Specific physics 
analyses will use selection criteria of higher efficiency and will achieve better angular resolutions 
than the ~ 2°, shown here. 
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Figure 11. The dependence of the space angle deviation of the LH fit on the muon energy for 
AMANDA-II. Shown are mean (stars) and median (circles) for simulated atmospheric neutrinos. 
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Figure 12. The space angle deviations of the LH fit as a function of the cosine of the incident 
zenith angle (for AMANDA-II). Shown are the mean (stars) and median (circles) for simulated 
atmospheric neutrinos. 

dependent shift of the reconstructed zenith angle and no systematic shift in azimuth. 
This is shown in figure 13 for simulated atmospheric neutrinos in AMANDA-II. The size 
of this shift depends on the zenith angle itself, and it is determined by the geometry 
of AMANDA, which has a larger size in vertical than in horizontal directions. From a 
comparison with AMANDA-BIO data [46,56], we observe that these shifts become smaller 
with a larger horizontal detector size. These shifts are confirmed by analyzing AMANDA 
events coincident with SPASE (see below). 

These angular deviations have been obtained from Monte Carlo simulations. They can be 
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Figure 13. The zenith angle shift of the reconstruction versus the cosine of the incident angle. 
Shown are mean (stars) and median (circles) for simulated atmospheric neutrinos. 




Space angle (degrees) 

Figure 14. Distribution of the space angle deviations between air shower directions assigned 
by SPASE-2 and muon directions assigned by AMANDA-BIO for coincident events measured in 
1997. The figure is not corrected for the systematic shift. 

experimentally verified by analyzing coincident events between AMANDA and SPASE. An 
analysis of data from the 10 string AMANDA-BIO detector, shown in figure 14, confirms 
the estimate of ~ 3° obtained from Monte Carlo studies for AMANDA-BIO . Unfolding 
the estimated SPASE resolution of ~ 1° confirms the estimated AMANDA-BIO resolution 
of ~ 3° near the SPASE-AMANDA coincidence direction [8-10]. 

A simulation-independent estimate can be obtained by splitting the hits of individual 
events in two parts and reconstructing each sub-event separately. The difference in the 
two results gives an estimate of the total angular resolution. Such analyses are being 
performed at present and results will be published separately. 
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7.3 Energy Reconstruction 
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Figure 15. Comparison of the resolution of three different energy reconstruction approaches for 
AMANDA-BIO. E gen is the generated energy (MC) and E iec the reconstructed energy. 

The energy resolution of the three methods, described in section 3.2.4, is shown in figure 15 
as function of the muon energy at its closest point to the AMANDA-BIO center. The 
resolution for AMANDA-BIO in Alog 10 i? is ~ 0.4, for the interesting energy range of a 
few TeV to 1 PeV. Below ~ 600 GeV the energy resolution is limited, because the amount 
of light emitted by a muon is only weakly dependent on its energy. Above 1 TeV the 
resolution improves because radiative energy losses become dominant. Above 100 TeV the 
resolution degrades, because energy loss fluctuations dominate. 

Although these methods are quite different, their performances are similar. The full E reco 
and Phu methods achieve similar resolutions up to 1 PeV. The Phu method becomes worse 
above this energy, because in AMANDA-B10 almost all of the OMs are hit, and the method 
saturates. In contrast, the Neural Net method shows a slightly poorer resolution up to 
1 PeV but is better above. Its resolution is relatively constant over several decades of 
energy. This is an advantage when reconstructing an original energy spectrum with an 
unfolding procedure as in [36]. 

The AMANDA-II detector contains more than twice as many OMs as AMANDA-B10, and 
the energy resolution is better, especially at larger energies, er(A log 10 E) ~ 0.3. The neural 
net reconstruction results for AMANDA-II are shown in figure 16. Finally, the recently in- 
stalled transient waveform recorders (TWR) allow better amplitude measurements, which 
should significantly improve the results of the energy reconstructions, in particular, the 
full E reco method [57]. 

As discussed in section 1, the cascade channel can achieve substantially better resolutions, 
because the full energy is deposited inside or close of the detector. Energy resolutions 
in Alog 10 £ of < 0.2 and < 0.15 can be achieved by AMANDA-B10 and AMANDA-II, 
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Figure 16. Energy reconstruction for simulated muons of different fixed energy in AMANDA-II, 
using the neural net method. 

respectively [4]. 



7.4 Systematic Uncertainties 
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Figure 17. The zenith angle deviations (RMS) as function of an additional uncertainty in the 
to time calibration. Data is shown for simulated atmospheric neutrino events in AMANDA-BIO 
with the selection of [30]. The transit time of the PMTs has been shifted without correcting for 
in the reconstruction. The shift is a fixed value for each PMT, obtained from a random Gaussian 
distribution. 

Several parameters of the detector are calibrated and therefore only known with limited 
accuracy. These parameters include the time offsets, the OM positions and the absolute 
OM sensitivities. We have estimated the effects of these uncertainties on the resolution of 
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AMANDA reconstructions [56]. As an example, figure 17 shows the effect of an additional 
contribution to the time calibration uncertainty for the 10 string AMANDA-BIO detector. 
The zenith angular resolutions for simulated atmospheric neutrino events only degrade 
when the additional timing uncertainties exceed 10 ns. Additional tests with similar re- 
sults were done with non-random systematic shifts such as a depth dependent shift or a 
string dependent shift. Therefore, the angular resolution is insensitive to the uncertain- 
ties in the time calibration. The geometry of the detector is known to better than 30 cm 
horizontally and to better than 1 m vertically, which corresponds to timing uncertain- 
ties of < 1 or 3.5 ns, respectively. Therefore, the geometry calibration is also sufficiently 
accurate. 

Similarly, the effect of uncertainties on other parameters, like the absolute PMT effi- 
ciency, has been investigated. No indication was found that the remaining calibration 
uncertainties seriously affect the angular resolution or the systematic zenith angle off- 
set. The combined calibration uncertainties are expected to affect the accuracy of the 
reconstruction by less than 5 % in the zenith angle resolution and to less than 0.5° in the 
absolute pointing offset. 



8 Discussion and Outlook 

We have developed methods to reconstruct and identify muons induced by neutrinos [30], 
inspite of the challenges of the natural environment and large backgrounds. These meth- 
ods allow us to establish AMANDA as a working neutrino telescope. The reconstruction 
techniques described in this paper are still subject to improvement in several aspects: 

The likelihood description: The likelihood functions for track reconstruction are based 
on the assumption of exactly one infinitely long muon track per event. Extensions of 
this model to encompass starting muon tracks (including the description of the hadronic 
vertex), stopping muons, muon bundles of non- negligible width, and multiple indepen- 
dent muons will be important, particularly in the context of larger detectors such as Ice 
Cube. Initial efforts fitting multiple muons with the direct walk algorithm have been 
useful in rejecting coincident down-going muons, and work toward reconstructing muon 
bundles has begun in the context of events coincident with SPASE air showers. 

The p.d.f. calculation: The likelihood function is based on parametrizations of prob- 
ability density functions (p.d.f.). The p.d.f. is obtained from Monte Carlo simulations, 
and its accuracy is limited by the accuracy of the simulation. Better simulations lead 
directly to a better p.d.f. and hence better reconstructions. 

The p.d.f. parametrizations: The p.d.f. is parametrized by functions (e.g. the Pandel 
functions) which only approximate the full p.d.f.. More accurate parametrization func- 
tions will result in better reconstructions. For example, the scattering coefficient shows 
a significant depth dependence (see section 2). The current reconstruction is based on 
an average p.d.f. assuming depth-independent ice properties. While the track recon- 
struction is relatively insensitive to the accuracy of the parametrization, we expect a 
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depth-dependent p.d.f. to have better energy reconstruction. 

Complementary information: The current reconstruction algorithms do not include 
all available information in an event. In particular, correlations between detected PMT 
signals are ignored. For this reason dedicated selection parameters have been designed to 
exploit this information. They are used to discriminate between well reconstructed and 
poorly reconstructed events and improve the quality of the data sample. Future work 
will try to improve these parameters and expand the present likelihood description. 

Transient waveform recorders: At the beginning of the year 2003, the detector read- 
out has been upgraded with transient waveform recorders [57]. We expect a substantial 
improvement of the multiple-photon detection and the dynamic range in particular for 
high muon energies. 

The construction of a much larger detector, the IceCube detector, will start in the year 
2004. It will consist of 4800 PMT deployed on 80 vertical strings and will surround the 
AMANDA detector [58]. The performance of IceCube has been studied with realistic Monte 
Carlo simulations and similar analysis techniques as described in this paper [59]. The 
result is a substantially improved performance in terms of sensitivity and reconstruction 
accuracy. A direction accuracy of about 0.7° (median) for energies above 1 TeV is achieved. 
Similar to AMANDA, we expect a further improvement by exploiting the full information, 
available from the recorded wave-forms, in the reconstruction. 
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